This project aimed to look at the spatial variability of color morph and Symbiodinium clades C and D present in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distributions at scales ranging from location within the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch) and depths. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of potential stress-response in Kane’ohe Bay.
knitr::opts_knit$set(root.dir = normalizePath(".."))
library(data.table)
library(effects)
library(lsmeans)
library(devtools)
library(plyr)
library(reshape2)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
library(png)
library(ecodist)
library(nnet)
# Import Sample Metadata
Coral_Data <- read.csv("Data/Collection Data/Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
# Import qPCR Data
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
## SHA-1 hash of file is 5d640f385ee8e972e6c6e54b60abb88689e04046
Mcap.plates <- list.files(path = "Data/qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm = list(C=2.26827, D=0),
copy.number=list(C=33, D=3))
Mcap <- Mcap$result
## Remove positive and negative controls from qPCR data
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
## Change "Sample.Name" column to "Colony"
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
## Filter out detections without two technical replicates
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
## Set C:D ratio to -Inf (when only D) and Inf (when only C)
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
## Order qPCR data by Colony
Mcap <- Mcap[with(Mcap, order(Colony)), ]
## Calculate Proportion C and Proportion D
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
## Categorize each sample as C-dominated or D-dominated
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
## Define category for symbiont mixture
Mcap$Mix <- factor(ifelse(Mcap$propC>Mcap$propD,
ifelse(Mcap$propD!=0, "CD", "C"),
ifelse(Mcap$propD>Mcap$propC,
ifelse(Mcap$propC!=0, "DC", "D"), NA)),
levels = c("C", "CD", "DC", "D"))
# Merge Sample Metadata with qPCR Data
Symcap <- merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
# Import Color Data
colscores <- read.csv("Data/Collection Data/Color_Scores.csv")
## Determine color scores by majority rule, and number of agreements
colscores$majority <- apply(colscores[,-1], 1, function(x) names(table(x))[which.max(table(x))])
colscores$n.agree <- apply(colscores[,-c(1,7)], 1, function(x) max(table(x)))
## Set which set of color scores will be used in analysis (=majority)
Symcap <- merge(Symcap, colscores[,c("Colony","majority","n.agree")], all=T)
Symcap$Color.Morph <- factor(Symcap$majority)
## Set color/symbiont groupings
Symcap$ColDom <- interaction(Symcap$Color.Morph, Symcap$Dom)
To account for the influence of tides, depth was adjusted according to the differences in mean sea level from the recorded sea level on each collection day to the nearest 6-minute interval. Mean sea level was obtained from NOAA daily tide tables for Moku o Lo’e.
JuneTide=read.csv("Data/Tide_Tables/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Data/Tide_Tables/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("Data/Tide_Tables/Station_1612480_tide_ht_20160801-20160812.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)),
format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2
# Add estimated times for missing values
Symcap$Time[170] <- as.POSIXct("2016-06-14 12:07:00")
Symcap$Time[177] <- as.POSIXct("2016-06-14 12:20:00")
Symcap$Time[178] <- as.POSIXct("2016-06-14 12:22:00")
Symcap$Time[180] <- as.POSIXct("2016-06-14 13:08:00")
Symcap$Time[187] <- as.POSIXct("2016-06-14 12:42:00")
Symcap$Time[188] <- as.POSIXct("2016-06-14 12:44:00")
Symcap$Time[206] <- as.POSIXct("2016-06-16 13:10:00")
Symcap$Time[208] <- as.POSIXct("2016-06-16 13:24:00")
Symcap$Time[211] <- as.POSIXct("2016-06-16 12:37:00")
Symcap$Time[218] <- as.POSIXct("2016-06-16 12:27:00")
Symcap$Time[448] <- as.POSIXct("2016-07-16 13:32:00")
Round6 <- function (times) {
x <- as.POSIXlt(times)
mins <- x$min
mult <- mins %/% 6
remain <- mins %% 6
if(remain > 3L) {
mult <- mult + 1
} else {
x$min <- 6 * mult
}
x <- as.POSIXct(x)
x
}
Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
Symcap<-merge(Symcap, Tide, by="Time.r", all.x=T)
Symcap$newDepth <- Symcap$Depth..m.- Symcap$TideHT
Symcap$DepthInt <- cut(Symcap$newDepth, breaks = 0:13)
Symcap.f <- Symcap[which(Symcap$n.agree > 3), ]
Collection sites represented on this map indicate the 9 fringe and 16 patch reef sites for a total of 25 collection reefs. In total, 707 colonies were collected for the analysis.
# Define map area and get satellite data from Google using RgoogleMaps
#KB <- c(21.46087401, -157.809907)
#KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = #FALSE)
#save(KBMap, file = "Data/KBMap.Rdata")
load("Data/KBMap.Rdata")
# Get representative coordinates for each reef based on sample GPS data
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
# Remove mislabeled reef causing duplicate colony ID's
XY <- subset(XY, Reef.ID!="37")
XY <- merge(XY, unique(Symcap[, c("Reef.ID", "Reef.Type", "Bay.Area")]), by="Reef.ID")
# Plot KBay Map
png(filename = "Figures/Fig1.png", width = 169, height = 120, units = "mm", res=300)
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude,
pch=c(1,2)[XY$Reef.Type],
col=c("red", "green", "blue")[XY$Bay.Area],
lwd=2, cex = 1.2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
legend("bottomleft", legend=c("Northern", "Central", "Southern", "Fringe", "Patch"),
pch=c(22, 22, 22, 1, 2), col=c("green", "red", "blue", "black", "black"),
pt.bg=c("green", "red", "blue", "white", "white"))
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("Data/coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()
dev.off()
Figure 1
Number of colonies dominated by either clade C or clade D.
table(Symcap$Dom)
##
## C D
## 434 273
Number of colonies with clade C only (C), C in greater abundance than D (CD), D in greater abundance than C (DC), or clade D only (D).
table(Symcap$Mix)
##
## C CD DC D
## 376 58 264 9
The proportion of clade C- and D-dominated corals with a single symbiont vs. mixtures.
results <- table(Symcap$Mix, Symcap$Dom)
out <- prop.table(results, margin = 2)
names(dimnames(out)) <- c("composition", "dominant symbiont")
out
## dominant symbiont
## composition C D
## C 0.86635945 0.00000000
## CD 0.13364055 0.00000000
## DC 0.00000000 0.96703297
## D 0.00000000 0.03296703
png("Figures/Fig3.png", width=81, height=81, units="mm", res=300)
par(mar=c(3, 3, 1, 4.5), mgp=c(2,0.5,0), tcl=-0.25, cex=0.8)
barplot(prop.table(results, margin = 2), col = c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), xlab = "Dominant Symbiont", ylab = "Proportion of Colonies")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), inset = c(-.35, 0), xpd = NA)
dev.off()
## quartz_off_screen
## 2
Figure 3
“D-only” hardly ever occurs. In fact, when looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.
The overall prevalence of orange and brown colonies in the dataset. Note that for all analyses involving color, colonies with amiguous color assignemnts (only 3/5 observers assigning the same color morph) are filtered out of the dataset.
table(Symcap.f$Color.Morph)
##
## Brown Orange
## 254 287
The number (first table) and proportion (second table) of brown and orange colonies that were dominated by clade C or D. A chi-squared test was used to test for differences in these proportions.
results <- table(Symcap.f$Dom, Symcap.f$Color.Morph)
results
##
## Brown Orange
## C 246 70
## D 8 216
prop.table(results, margin = 2)
##
## Brown Orange
## C 0.96850394 0.24475524
## D 0.03149606 0.75524476
chisq.test(results)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 287.32, df = 1, p-value < 2.2e-16
propD <- Symcap.f$propD[which(Symcap.f$propD > 0 & Symcap.f$propD < 1)]
propDHist <- subset(Symcap.f, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(Symcap.f, propD==0), table(Color.Morph))
o <- with(subset(Symcap.f, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
png("Figures/Fig4.png", width=81, height=81, units="mm", res=300)
par(mar=c(3, 3, 1, 0), cex=0.8, mgp=c(2,0.5,0), tcl=-0.25, cex.axis=0.7)
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0%", "20%", "40%", "60%", "80%", "<100%"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)
dev.off()
## quartz_off_screen
## 2
Figure 4
A logistic regression is used to model the probability of a colony being D-dominated as a function of depth, bay area (N, S, Central), and reef type (patch, fringe).
# Fit GLM for dominant symbiont
Symcap$Dominant <- ifelse(Symcap$Dom=="C", 0, 1)
dom <- glm(Dominant ~ newDepth * Bay.Area * Reef.Type, family="binomial", Symcap)
anova(dom, test="Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 705 942.15
## newDepth 1 129.043 704 813.10 < 2.2e-16 ***
## Bay.Area 2 1.268 702 811.84 0.530554
## Reef.Type 1 0.050 701 811.79 0.823100
## newDepth:Bay.Area 2 0.874 699 810.91 0.646002
## newDepth:Reef.Type 1 8.005 698 802.91 0.004665 **
## Bay.Area:Reef.Type 2 5.142 696 797.76 0.076442 .
## newDepth:Bay.Area:Reef.Type 2 1.649 694 796.12 0.438459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The model shows that depth is the most significant predictor, but that the effect of depth is different for patch vs. fringe reefs, due to the significant Depth * Reef.Type interaction. We plot the overall effect of depth averaged across the entire bay, as well as the specific effect of depth for patch reefs and fringing reefs in the figure below.
# Get fitted values for overall depth effect
dom.fit <- data.frame(Effect(focal.predictors = "newDepth", mod = dom,
xlevels = list(newDepth=seq(0,11,0.1))))
dom.fit2 <- data.frame(Effect(focal.predictors=c("newDepth", "Reef.Type"), mod = dom,
xlevels = list(newDepth=seq(0,11,0.1))))
# Break depth into intervals for plotting
Symcap$DepthInt <- cut(Symcap$newDepth, breaks = 0:13)
Symcap$Dominant2 <- ifelse(Symcap$Dom=="C", 1, 0)
results=table(Symcap$Dominant2, Symcap$DepthInt)
dom.props <- prop.table(results, margin = 2)
# Plot
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(dom.props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D", "KBay", "patch", "fringe"),
pch=c(22, 22, NA, NA, NA),
pt.bg=c(alpha("blue", 0.25), alpha("red", 0.25), NA, NA, NA),
xpd = NA, lty=c(0,0,1,3,2), pt.cex=c(2,2,1,1,1), lwd=c(NA,NA,2,1,1))
par(new = T)
par(mar=c(4.2, 4, 2, 6))
plot(dom.fit$fit~dom.fit$newDepth, xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of D-dominance")
# Depth at which probability of C vs. D dominance is 50/50
seq(0,11,0.1)[which.min(abs(dom.fit$fit-0.5))]
## [1] 1.3
with(subset(dom.fit2, Reef.Type=="Fringe"), lines(fit ~ newDepth, lty=2))
with(subset(dom.fit2, Reef.Type=="Patch"), lines(fit ~ newDepth, lty=3))
Bars indicate relative frequency of occurence of clade C vs. D dominance at 1m depth intervals when pooling all samples collected. Overlaid lines indicate probability of a colony being D-dominated for the entire bay (thick solid line), and for fringing reefs (dashed line) and patch reefs (dotted line).
# Fit GLM for color morph
Symcap.f$Color <- ifelse(Symcap.f$Color.Morph=="Orange", 0, 1)
cm <- glm(Color.Morph ~ newDepth * Bay.Area * Reef.Type, family="binomial", Symcap.f)
anova(cm, test="Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Color.Morph
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 540 747.97
## newDepth 1 65.959 539 682.01 4.604e-16 ***
## Bay.Area 2 1.281 537 680.73 0.52709
## Reef.Type 1 0.735 536 680.00 0.39112
## newDepth:Bay.Area 2 1.622 534 678.37 0.44440
## newDepth:Reef.Type 1 4.859 533 673.51 0.02750 *
## Bay.Area:Reef.Type 2 1.210 531 672.31 0.54616
## newDepth:Bay.Area:Reef.Type 2 4.792 529 667.51 0.09108 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Get fitted values
cm.fit <- data.frame(Effect(focal.predictors = "newDepth", mod = cm,
xlevels = list(newDepth=seq(0,11,0.1))))
cm.fit2 <- data.frame(Effect(focal.predictors =c("newDepth", "Reef.Type"), mod = cm,
xlevels = list(newDepth=seq(0,11,0.1))))
# Plot
results <- table(Symcap.f$Color, Symcap.f$DepthInt)
cm.props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(cm.props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange", "KBay", "patch", "fringe"),
pch=c(22, 22, NA, NA, NA),
pt.bg=c(alpha("sienna", 0.25), alpha("orange", 0.25), NA, NA, NA),
xpd = NA, lty=c(0,0,1,3,2), pt.cex=c(2,2,1,1,1), lwd=c(NA,NA,2,1,1))
par(new = T)
par(mar=c(4.2, 4, 2, 6))
plot(cm.fit$fit~cm.fit$newDepth, xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of Orange-Dominance")
# Depth at which probability of orange or brown is 50/50
seq(0,11,0.1)[which.min(abs(cm.fit$fit-0.5))]
## [1] 2.7
with(subset(cm.fit2, Reef.Type=="Fringe"), lines(fit ~ newDepth, lty=2))
with(subset(cm.fit2, Reef.Type=="Patch"), lines(fit ~ newDepth, lty=3))
Bars indicate relative frequency of brown vs. orange color morph dominance at 1m depth intervals when pooling all samples collected. The overlaid line indicates the probablity of a colony to be orange across depths.
While brown colonies were always more likely C-dominated, orange colonies were observed to switch from D to C dominance. The depth threshold where this switch occurred is represented here.
# Probability of D dominance in Orange and Brown colonies across Depth
Symcap.f$Dominant <- ifelse(Symcap.f$Dom=="C", 0, 1)
mod <- glm(Dominant ~ newDepth * Color.Morph, family="binomial", data=Symcap.f)
anova(mod, test="Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 539 732.85
## newDepth 1 99.503 538 633.35 < 2e-16 ***
## Color.Morph 1 277.659 537 355.69 < 2e-16 ***
## newDepth:Color.Morph 1 6.219 536 349.47 0.01264 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test to see if relationship between depth and D-dominance in orange colonies differs between patch and fringe reefs (not enough Brown-D colonies to support full factorial model)
mod2 <- glm(Dominant ~ newDepth * Reef.Type, family="binomial",
data=subset(Symcap.f, Color.Morph=="Orange"))
anova(mod2, test="Chisq") # depth*reef.type interaction is NOT significant
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 285 318.32
## newDepth 1 39.879 284 278.44 2.701e-10 ***
## Reef.Type 1 0.047 283 278.39 0.82770
## newDepth:Reef.Type 1 3.525 282 274.87 0.06047 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Get fitted values for probability of D-dominance in orange and brown colonies across depth
int.fit <- data.frame(Effect(focal.predictors=c("newDepth", "Color.Morph"), mod=mod,
xlevels=list(newDepth=seq(0,11,0.1)), type="response"))
with(subset(int.fit, Color.Morph=="Orange"), {
plot(fit ~ newDepth, ylim = c(0,1), type="l", col="orange",
lwd=3, xlab="", ylab="")
})
with(subset(int.fit, Color.Morph=="Brown"), {
lines(fit ~ newDepth, col="sienna", lwd=3)
})
abline(h = 0.5, lty=2)
# Point at which probability of C or D dominance is 50/50
with(subset(int.fit, Color.Morph=="Orange"), seq(0,11,0.1)[which.min(abs(fit-0.5))])
## [1] 4.3
# (not applicable for Brown colonies)
# Probability of D-dominance at 1m
with(subset(int.fit, Color.Morph=="Orange"), fit[seq(0,11,0.1)==1])
## [1] 0.8377806
with(subset(int.fit, Color.Morph=="Brown"), fit[seq(0,11,0.1)==1])
## [1] 0.03358084
# Probability of D-dominance at 10m
with(subset(int.fit, Color.Morph=="Orange"), fit[seq(0,11,0.1)==11])
## [1] 0.03694416
with(subset(int.fit, Color.Morph=="Brown"), fit[seq(0,11,0.1)==11])
## [1] 0.02537863
This three-panel figure summarizes the influence of depth (and interaction with reef type) on the dominant symbiont (top) and color morph of M. capitata (middle), and the differing influence of depth on dominant symbiont in brown vs. orange colonies (bottom).
## quartz_off_screen
## 2
Figure 5
A multinomial logistic regression is used to model the occurrence of all four color-symbiont combinations across depth for all samples from Kaneohe Bay.
prop.table(table(Symcap.f$ColDom))
##
## Brown.C Orange.C Brown.D Orange.D
## 0.45555556 0.12962963 0.01481481 0.40000000
depthmod <- multinom(ColDom ~ newDepth, Symcap.f)
## # weights: 12 (6 variable)
## initial value 748.598955
## iter 10 value 533.941343
## iter 20 value 513.655715
## final value 513.374759
## converged
depthdat <- data.frame(newDepth = seq(0,12,0.5))
depthpred <- predict(depthmod, newdata=depthdat, "probs")
stackpoly(depthpred[,c(4,2,3,1)], stack=T, col=rev(c("burlywood4","burlywood3","darkorange3","darkorange")),
ylim=c(0,1), xaxlab=depthdat$newDepth, xlab="Depth", ylab="Probability")
legend("topright", legend=c("Brown.C", "Brown.D", "Orange.C","Orange.D"),
fill=c("burlywood4","burlywood3","darkorange3","darkorange"))
A multinomial logistic regression is used to model the occurrence of all four color-symbiont combinations across depth for patch reefs and fringing reefs separately (based on significant interaction between depth and reef type).
## Fit multinomial model
mod <- multinom(ColDom ~ newDepth * Bay.Area * Reef.Type, Symcap.f)
## # weights: 52 (36 variable)
## initial value 748.598955
## iter 10 value 520.226537
## iter 20 value 496.986515
## iter 30 value 492.684338
## iter 40 value 492.166498
## iter 50 value 492.136316
## final value 492.136065
## converged
mod <- step(mod, trace=0)
## trying - newDepth:Bay.Area:Reef.Type
## # weights: 44 (30 variable)
## initial value 748.598955
## iter 10 value 526.930615
## iter 20 value 498.413129
## iter 30 value 494.709358
## iter 40 value 494.514392
## iter 50 value 494.507500
## final value 494.506893
## converged
## trying - newDepth:Bay.Area
## trying - newDepth:Reef.Type
## trying - Bay.Area:Reef.Type
## # weights: 36 (24 variable)
## initial value 748.598955
## iter 10 value 521.524472
## iter 20 value 499.504913
## iter 30 value 498.459588
## iter 40 value 498.431873
## final value 498.430819
## converged
## trying - newDepth:Reef.Type
## trying - Bay.Area:Reef.Type
## # weights: 28 (18 variable)
## initial value 748.598955
## iter 10 value 524.918365
## iter 20 value 502.712790
## iter 30 value 502.002427
## final value 502.002422
## converged
## trying - Bay.Area
## trying - newDepth:Reef.Type
## # weights: 20 (12 variable)
## initial value 748.598955
## iter 10 value 511.828123
## iter 20 value 504.944395
## final value 504.776669
## converged
## trying - newDepth:Reef.Type
mod # Only depth:Reef.Type interaction is significant
## Call:
## multinom(formula = ColDom ~ newDepth + Reef.Type + newDepth:Reef.Type,
## data = Symcap.f)
##
## Coefficients:
## (Intercept) newDepth Reef.TypePatch newDepth:Reef.TypePatch
## Orange.C -1.297780 -0.06478447 0.351041 -0.005528741
## Brown.D -1.266232 -3.47581513 -2.485392 3.576803842
## Orange.D 1.459716 -0.88420417 -0.597042 0.452442790
##
## Residual Deviance: 1009.553
## AIC: 1033.553
## Get fitted values
depthdat <- expand.grid(newDepth=seq(0,10,0.5), Reef.Type=c("Patch", "Fringe"))
depthpred <- predict(mod, newdata=depthdat, "probs")
### Plotting
par(mfrow=c(2,1), mar=c(2,3,2,2), mgp=c(1.5,0.3,0), tcl=-0.25, cex.axis=0.8)
stackpoly(depthpred[1:21,c(4,2,3,1)], stack=T,
col=rev(c("burlywood4","burlywood3","darkorange3","darkorange")),
ylim=c(0,1), xaxlab=seq(0,10,0.5), xlab="", ylab="Probability")
title("A. Patch reefs", adj=0)
legend("topright", legend=c("Brown-C", "Brown-D", "Orange-C", "Orange-D"),
fill = c("burlywood4","burlywood3","darkorange3","darkorange"),
inset=c(0.2,0.05), cex=0.8, bty="n")
par(mar=c(3,3,2,2))
stackpoly(depthpred[22:42,c(4,2,3,1)], stack=T,
col=rev(c("burlywood4","burlywood3","darkorange3","darkorange")),
ylim=c(0,1), xaxlab=seq(0,10,0.5), xlab="Depth (m)", ylab="Probability")
title("B. Fringe reefs", adj=0)
Bray-Curtis dissimilarity among reefs is calculated based on the relative abundance of each color-symbiont phenotype at each reef (after correcting for the influence of depth across the whole bay to correct for differences in sampling effort across depth at different reefs). A Mantel test is then used to test whether Bray-Curtis dissimilarity is correlated with physical (Euclidean) distance, to test the alternative hypothesis that reefs closer together in space show more similar occurrences of color and dominant symbiont.
model <- multinom(ColDom~newDepth + Reef.ID, Symcap.f)
## # weights: 108 (78 variable)
## initial value 748.598955
## iter 10 value 489.030699
## iter 20 value 468.198426
## iter 30 value 467.753739
## iter 40 value 467.641730
## iter 50 value 467.550924
## iter 60 value 467.520408
## final value 467.519600
## converged
newdat <- data.frame(Reef.ID = levels(Symcap.f$Reef.ID),
newDepth = mean(Symcap.f$newDepth))
pred <- round(predict(model, newdata = newdat, "probs"), 7)
adjprops <- data.frame(Reef.ID=as.character(newdat[,1]), pred)
actprops <- data.frame(
Reef.ID=as.character(newdat[,1]),
as.data.frame.matrix(prop.table(table(Symcap.f$Reef.ID, Symcap.f$ColDom), margin=1)))
# Test for spatial autocorrelation in color-symbiont distributions
# Get Lat/Lon positions for each reef
Latitude=aggregate(Latitude~Reef.ID, data=Symcap.f, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap.f, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
# Mantel test comparing Bray-Curtis dissimilarity among reefs (composition of color-symbiont groups) vs. euclidean/physical distance among reefs
df <- merge(XY, adjprops, by="Reef.ID")
reef.dists <- dist(cbind(df$Longitude, df$Latitude))
bc.dists <- vegan::vegdist(df[,4:7], method="bray")
set.seed(12456)
ecodist::mantel(bc.dists~reef.dists) # No spatial autocorrelation
## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## 0.05647236 0.22100000 0.78000000 0.50000000 -0.04468309 0.13257895
No significant spatial auto-correlation is detected.
The depth-corrected relative abundances of each color-symbiont combination at each reef are plotted on a map of Kaneohe Bay.
png(filename = "Figures/Fig6.png", width = 169, height=120, units="mm", res=300)
par(oma=c(3,3,0,0), mfrow=c(1,1))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, pch=NA)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
newcoords <- LatLon2XY.centered(KBMap, df$Latitude, df$Longitude, zoom=13)
df$X <- newcoords$newX
df$Y <- newcoords$newY
rownames(df) <- df$Reef.ID
df <- df[, -1]
df <- na.omit(df)
df[df==0] <- 0.0001
apply(df, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=as.numeric(reef[c("Brown.C", "Orange.C", "Brown.D", "Orange.D")]),
radius = 8, col = c("burlywood4","darkorange3","burlywood3","darkorange"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("burlywood4","darkorange3","burlywood3","darkorange"))
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("Data/coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()
dev.off()
Figure 6
Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.